Transmission model

Author

Elena Colman & Fiona Harris

Transmission rate

What happens if the transmission rate is dependent on the age groups?

Model diagram: age-dependent transmission rates

To examine the impact of age-specific transmission dynamics, we modified the model to account for different contact patterns between juveniles and adults. Instead of using a single uniform transmission rate across all age groups, we introduced three distinct contact rates to reflect behavioral differences with following transmission rate values:

  • Juvenile-to-juvenile transmission (e.g., through play and nesting interactions): \(r_{\beta 1}\)
  • Adult-to-adult transmission (e.g., through breeding and fighting): \({r_\beta 2}\)
  • Cross-age transmission between juveniles and adults (less frequent interaction): \(r_{\beta 3}\)

These values allow us to explore how increased transmission among adults (e.g., due to aggressive behaviors) might influence infection dynamics in the overall population, and how juvenile–adult interactions contribute to disease spread.

Show code 
TransmissionModel <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    NJ <- SJ + EJ + IJ
    NAd <- SAd + EAd + IAd
    Y <- 1 / (1 + exp(-rho * (NJ + NAd - kappa)))
    X <- 1 - Y
    delta <- 0.1
    total_population <- NJ + NAd

    births_total <- X * lambda * (SAd + EAd + IAd)
    birth_SJ <- X * lambda * (SAd + EAd) + X * lambda * IAd * (1 - delta)
    births_EJ <- X * lambda * delta * IAd

    infectionJuveniles <- r_beta1 * (IJ / NJ)
    infectionAdults <- r_beta2 * (IAd / NAd)
    infectionCross <- r_beta3 * (IJ / NJ + IAd / NAd)

    dSJ <- birth_SJ - (mu * Y + gammma) * SJ - (infectionJuveniles + infectionCross) * SJ
    dEJ <- (infectionJuveniles + infectionCross) * SJ + births_EJ - (mu * Y + gammma + alpha) * EJ
    dIJ <- alpha * EJ - (mu * Y + gammma + tau) * IJ

    dSAd <- gammma * SJ - mu * Y * SAd - (infectionAdults + infectionCross) * SAd
    dEAd <- gammma * EJ + (infectionAdults + infectionCross) * SAd - (mu * Y + alpha) * EAd
    dIAd <- gammma * IJ + alpha * EAd - (mu * Y + tau) * IAd

    
      # Calculate deaths
    natural_deaths <- mu * Y * (SJ + EJ + IJ + SAd + EAd + IAd)
    tb_deaths <- tau * (IJ + IAd)
    
    # Return deaths as extra outputs
    return(list(c(dSJ, dEJ, dIJ, dSAd, dEAd, dIAd),
                natural_deaths = natural_deaths,
                tb_deaths = tb_deaths, total_population = total_population))}
    
)
}


# Simulation grid
yini <- c(SJ = 20, EJ = 0, IJ = 0, SAd = 30, EAd = 0, IAd = 1)
times <- seq(0, 7.5, by = 0.1)
params_fixed <- c(lambda = 5, kappa = 50, rho = 0.5, mu = 1/5, tau = 1, gammma = 1, alpha = 5)

# Generate simulation grid
param_grid <- expand_grid(
  r_beta1 = seq(0.5, 4.5, 0.5),
  r_beta2 = seq(0.5, 4.5, 0.5),
  r_beta3 = seq(0.5, 4.5, 0.5)
)

# Run model for each combo
sim_df <- param_grid |> 
  rowwise() |>
  mutate(sim = list(as.data.frame(
  ode(y = yini, times = times, func = TransmissionModel,
      parms = c(params_fixed, r_beta1 = r_beta1, r_beta2 = r_beta2, r_beta3 = r_beta3))
))) |>
  unnest(sim) |>
  ungroup()

ojs_define(possum_data = purrr::transpose(as.data.frame(sim_df)))


maxPopTransmission <- max(sim_df$total_population)
timeAtMaxPopTransmission <- sim_df$time[which.max(sim_df$total_population)]

ojs_define(
  maxPopTransmission = round(maxPopTransmission),
  timeAtMaxPopTransmission = round( timeAtMaxPopTransmission,1))